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Abstract 

A theoretical description of shell structure for charged particles in a harmonic trap is explored 
at strong coupling conditions of T = 50 and 100. The theory is based on an extension of the 
hypernetted chain approximation to confined systems plus a phenomenological representation of 
associated bridge functions. Predictions are compared to corresponding Monte Carlo simulations 
and quantitative agreement for the radial density profile is obtained. 
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Systems of harmonically trapped charged particles exhibit a shell structure in their radial 
density profile. Recent studies [1, , both experimental and molecular dynamics 



simulation, indicate the localization of the charges on the surfaces of concentric spheres 
with a crystal-like ordering on each surface depending on the particle number N. These are 
effectively the zero temperature ground state for the system or energetically close excited 
states. An accurate and instructive theoretical description of these spherical crystals has 
been given recently using shell models 0, fl g A complementary theoretical description 
of the fluid phase at finite temperatures has now been given as well |9[, showing a closely 
related "blurred" shell structure that sharpens at stronger Coulomb coupling T = q 2 /r^ksT 
(q is the charge, r is the ion sphere radius in terms of the average trap density). This fluid 
phase theory demonstrates that correlations play the dominant role in the formation of shell 
structure. An extension of the hypernetted chain approximation (HNC) for bulk fluids to 
the localized charges of a trap was shown to provide all of the qualitative features of shell 
structure (e.g., number and location of density peaks, shell populations) for the Coulomb 
coupling constant T > 10 . However, comparison with Monte Carlo (MC) simulation results 
at selected conditions showed significant errors in the HNC amplitudes and widths of the 
shells (20 — 40%). An adjusted HNC (AHNC) was proposed to correct these deficiencies, 
using a phenomenological representation of the neglected "bridge functions". Excellent 
agreement with MC results was obtained in this way for T = 20, N = 25, 100 and T = 
40, N = 300. The objective here is to report further exploration of the AHNC at the 
stronger coupling values T = 50 and 100, and to note some interesting similarities to the 
crystal shell structure. 

The origin of the AHNC theory for confined systems is summarized briefly first. Consider 
a system of iV classical charges in an external potential. The Hamiltonian for this system is 

N N 1 1 N 

H = H + Y, Vo^), H = Y; ^v 2 + g E V ( r v) W 

i=l i=l i^ij=l 

where r« and are the position and velocity of charge i. The external potential seen by each 
particle is denoted by Vb(r), and the interaction between the pair i,j is V(rij) (application 
here will be limited to Coulomb interactions but the discussion is more general) . The external 
potential induces a non-uniform equilibrium density n(r). It follows from density functional 
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theory that n(r) obeys the equation [1 

^ = ^ (I) . ? «WM, (2) 

z on (r) 

where z = e^, /x is the chemical potential, and A = (h 2 f3/27rm) 1 ^ 2 is the thermal wavelength. 
The excess free energy F ex ((3 \ n) is a universal functional of the density for the Hamiltonian 
ifoj independent of the applied external potential Vo, and describes all correlations among 
the particles. The solutions to d2J) are such that there is a unique equilibrium density n(r) 
for each Vo(r), using the same F ex (/3 | n). A special case is the uniform density of a one 
component plasma (OCP), resulting from the Vo (r) for a uniform neutralizing background. 
The quite different non-uniform density of interest here results from the harmonic trap 
Vo (r) = muj 2 r 2 /2. (The spherical symmetry of Vo in both cases yields a spherically symmet- 
ric density profile in the fluid phase; the broken symmetry crystal profile is not considered 
here) . 

This observation that the OCP and trap densities are determined from the same excess 
free energy functional opens the possibility of describing correlations for the trap in terms 



of those of the OCP. This is done in reference 



9| with the result 



n(r)X 
In 



— [3-mu 2 r 2 + (3 I dr n(r')cocp (r — r\; T) + 5x(r|n) (3) 

Z A J 

where cocp is the direct correlation function of a one-component plasma, and B^{r\n) is the 
"bridge function" for the trapping potential. The Ornstein-Zernicke equation relates cocp 
to the pair distribution function for the OCP, <?ocp- A similar analysis gives an equation for 
9ogp 



In socp(r) = -Pq 2 ^ 1 + (3 J d? [<?ocp(r') - 1] c OC p (\r-^\;T) + £ocp(r|<7ocp)- (4) 

These two equations provide a formally exact description of the charged particle system 
from which approximations can be made. In particular, the HNC approximation for both 
the OCP and the trap density is defined by the neglect of both bridge functions in these 
equations. The trap density is then determined entirely in terms of correlations for the OCP. 

The HNC density profile is very accurate at weak to moderate coupling, F < 10, but 
has significant quantitative errors as the shell structure develops at larger values of the 
coupling. This is illustrated in Figured] for T = 30, N = 100. It was proposed in reference 
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FIG. 1: Comparison of HNC and adjusted HNC (AHNC) to Monte Carlo (MC) results. Conditions 
are r = 30 and N = 100. 

9j to improve the HNC by a phenomenological representation of the bridge functions in the 
form 

B(r\n) - \(T)(3V (r), (5) 

where Vq is the harmonic potential for or Coulomb potential for -Bocp- Also A(T) is 
chosen to interpolate between A(0) = and some constant A(oo). This approach was first 



introduced by Ng |ljj for the OCP pair distribution. Subsequently, Rosenfeld and Ashcroft 
have defined related "modified" HNC approximations using the bridge function as a fitting 
function [12j]. The advantage of the form ([5]) is that the effect of the bridge functions is to 
simply " renormalize" the external potential in both eqs. ([3]) and (j3J), so that the original 
HNC approximation is regained except with an effective coupling constant V defined by 

r' = [i + A(r)]r (6) 

In the original work of Ng, he obtained agreement with the Monte Carlo data for gocp to 
within a few percent using A(oo) = 0.6. That same value for A(oo) has been used in reference 
9| and in the results presented here. 

The improvement gained from this AHNC is also shown in Figure 1. All three curves agree 




FIG. 2: Comparison of density profile for AHNC and MC, for r = 50, N = 50 and 400 

regarding the number of peaks and the location of the peaks. Also, there is quantitative 
agreement regarding the number of charges in each shell (not shown). However, the HNC and 
AHNC results give very different heights and widths of the peaks themselves. Comparison 
to the MC results makes it clear that the bridge function effects are important at this value 
of the coupling constant. Similar improvement has been illustrated for T = 20 and 40 for 
several values of A" [9|. 

It remains to explore how well this approach extends to still larger values of T, motivated 
by its possible utility to study freezing into a spherical crystal. As a first step in that 
direction results are presented and discussed here for T = 50, 100. Figure ([2]) shows the 
density profile at T = 50 for two quite different particle numbers A^ = 50 and A" = 400. For 
small enough particle number, A" < 20, there is only one shell of finite radius near the point 
at which the density vanishes jjj. The latter is fixed by the condition that the Coulomb 
force on a test charge balances that of the trap, i.e. r max = N 1 ^ 3 . A second shell is present 
sA N — 50, while three additional shells are observed for A" = 400 (as well as a population 
at the origin [4]). The agreement between MC and AHNC is seen to be quite good in both 
cases except at very small r for which the both HNC and AHNC fail even at moderate 
coupling. 
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FIG. 3: Comparison of density profile for AHNC and MC results for T = 100 and N = 100. 

The discrepancies in amplitudes increase somewhat at very strong coupling. This is 
illustrated in Figure ([3]) showing the comparison of MC and AHNC results for T = 100, N = 
100. The error in peak height grows from about 10% for the outer shell to about 30% for 
the inner shell. Still the overall agreement remains quite good. To explore the conditions 
for freezing into spherical crystals, even stronger coupling conditions, T ~ 200 are expected 
to be relevant. However, several important features of the crystal ground state are already 
evident from the fluid phase AHNC results at the coupling conditions studied. Figure (jl]) 
shows a comparison of shell populations from ground state annealed MD simulations 2] with 
those from AHNC. The latter have been computed for 20 < T < 100, showing a very weak 
dependence on coupling strength. This is somewhat unexpected since the amplitudes and 
widths of the shells are quite sensitive to T. The shell populations and appearance of shells 
is clearly quite similar in the fluid and crystal phases. Figure ([5]) shows a corresponding 
comparison of shell locations as a function of N. Again the results are remarkably close for 
both fluid and crystal phases. 

The results for the ground state shown on Figures (jl]) and §5§ can be understood theoret- 
ically on the basis of a shell model. This represents the density as a sum of delta functions 
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FIG. 4: Comparison of shell populations from AHNC and MD as a function of particle number N. 
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The values of Thy , Tp cir6 determined by minimizing the total energy [6|, LZl] - A very accurate 
representation of the energy for this purpose is constructed on the basis of solutions to the 
Thomson problem Jjj] for ground state configurations of Coulomb charges on a single sphere 
8J. The intrashell energy is chosen as the known Thomson energy for each n u ,r u , while the 
intershell energy is taken as the monopole interactions for charges at the corresponding 
Thomson positions for each shell. In this way the crystal energy, shell populations, and 
locations are given very accurately. This opens the possibility to test free energies associated 
with ([3]) and (jl]) for broken rotational invariance using parameterized densities similar to 
the shell model (e.g., Gaussians centered at the Thomson positions) to identify a freezing 
transition. 

This work is supported by the Deutsche Forschungsgemeinschaft via SFB-TR 24, and by 
the NSF/DOE Partnership in Basic Plasma Science and Engineering under the Department 
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FIG. 5: Comparison of shell radii from AHNC and MD as a function of particle number N. 
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